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ABSTRACT 

In  this  paper  a  wide  class  of  difference  equations  is  described  for 
approximating  discontinuous  time  dependent  solutions,  with  prescribed 
initial  data,  of  hyperbolic  systems  of  nonlinear  conservation  laws.  Among 
these  schemes  we  determine  the  best  ones,  i.e.,  those  which  have  the  small¬ 
est  truncation  error  and  in  which  the  discontinuities  are  confined  to  a 
narrow  band  of  2-3  meshpoints.  These  schemes  are  tested  for  stability 
and  are  found  to  be  stable  under  a  mild  strengthening  of  the  Courant- 
%  Friedrichs -Lewy  criterion.  Test  calculations  of  one -dimensional  flows 

of  compressible  fluids  with  shocks,  rarefaction  waves  and  contact  discon¬ 
tinuities  show  excellent  agreement  with  exact  solutions.  In  particular, 
when  Lagrange  coordinates  are  used,  there  is  no  smearing  of  interfaces. 

The  additional  terms  introduced  into  the  difference  scheme  for  the 
purpose  of  keeping  the  shock  transition  narrow  are  similar  to,  although 
not  identical  with,  the  artificial  viscosity  terms,  and  the  like  of  them 
introduced  by  Richtmyer  and  von  Neumann  and  elaborated  by  other  workers 
in  this  field. 


* 
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1.  DIFFERENCE  SCHEMES  FOR  CONSERVATION  LAWS 

In  this  paper  we  consider  systems  of  conservation  laws,  i.e., 
equations  of  the  form 


\  -  V  d-D 

where  u  is  an  unknown  vector  function  of  x  and  t  with  n  components, 
and  f  a  given  vector  function  of  u,  depending  in  general  non-linearly 
on  u.  When  the  differentiation  on  the  right  side  of  (1.1)  is  carried  out 
a  quasi -linear  system  results: 

ut  =  Aux,  (1.1* ) 

where  A  =  A(u)  is  the  gradient  of  f .  We  require  that  our  system  of 
equations  he  hyperbolic  in  the  sense  that  the  matrix  A  has  n  distinct 
real  eigenvalues  p^,  . ..,  p^,  for  all  values  of  u.  Of  course  the 

eigenvalues  themselves  are  functions  of  u. 

The  negatives  of  the  quantities  p  are  the  local  sound  speeds. 

We  shall  index  them  in,  say,  monotonic  increasing  order. 
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The  initial  value  problem  is  to  find  a  solution  of  (1.1)  with 
prescribed  values  at  t  =  0: 

u(x,0)  =  cp(x).  (1.2) 

It  is  well  known  that  on  account  of  the  nonlinearity  of  (1.1)  no  smooth 
solution  will  in  general  exist  for  all  time.  Instead  we  have  to  seek 
weak  solutions  of  (1.1),  (1.2),  defined  by  the  requirement  that  the 
integral  relation 

(w  u  -  w  f )  dxdt  + 

"t  X 

be  satisfied  for  all  smooth  test  vectors  w  which  vanish  for  |x|  +  t 
large  enough. 

An  immediate  consequence  of  the  integral  relations  is  that  every 
piecewise  continuous  weak  solution  must  satisfy  the  Rankine-Hugoniot 
relation  across  a  line  of  discontinuity: 

s[u]  +  [f]  =  0;  (lA) 

here  the  brackets  denote  the  jump  across  the  discontinuity,  s  the  speed 
of  propagation  of  the  discontinuity. 

It  is  well  known  in  the  theory  of  systems  of  conservation  laws 
(see,  e.g.,  section  7  of  [9])  that  in  addition  one  has  to  impose  a  so- 


/ 


w(x,0)  cp(x)  dx  =  0 


(1.3) 
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called  entropy  condition  on  discontinuities.  This  condition  can  be 
formulated  as  follows: 

There  is  an  index  k  such  that  the  shock  speed  lies  between  the 
(k  -  1)  and  kth  characteristic  speeds  with  respect  to  the  state  on  the 
left  of  the  shock,  and  between  the  kth  and  (k  +  1)  characteristic 
speeds  on  the  right. 

We  shall  describe  now  a  fairly  large  class  of  difference  schemes 
that  can  be  used  to  approximate  weak  solutions  with  prescribed  initial 
data . 

First  we  choose  a  vector  valued  function  g  of  21  vector  argu¬ 
ments,  related  to  f  by  the  sole  requirement  that  when  all  the  21 
arguments  are  equal ,  g  reduces  to  f : 

g (u,  ...,  u)  =  f(u).  (1.5) 

Abbreviate  the  values  of  u(x)  at  the  lattice  points  x  +  kAx  and  at 
time  t  by  u^,  k  =  -  oo,  ...  oo.  We  define 

g(x+  Ax/2)  =  g(u_i+1,  u_i+2,  Uj); 

therefore  similarly  we  have 

g (x  -  Ax/2)  =  g(u_f,  u_f+1,  ...,  u^). 


« 
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We  take  now  the  following  difference  analogue  of  (1.1): 


Au  Ag 
At  =  Ax  ’ 


where  Au  denotes  the  forward  time  difference 


Au  =  u(x,  t+At)  -  u(x,t) 


and  Ag  the  symmetric  space  difference 


Ag  =  g(x  +  Ax/2)  -  g(x  -  Ax/2). 


(1.6) 


From  (1.6)  we  can  determine  u(x,  t+At): 

u(x,  t+At)  =  u(x,t)  +  XAg,  (1-6*) 

where  X  abbreviates  the  quotient  At/Ax  of  the  time  and  space  incre¬ 
ments.  Given  the  initial  values  cp  of  u,  we  can,  using  (1.6'), 
determine  successively  the  values  of  u  at  all  times  which  are  integer 
multiples  of  At.  We  claim  now  that  as  consequence  of  (1-5)  our  differ¬ 
ence  scheme  is  consistent  with  the  differential  equation  in  the  following 
sense:  denote  by  v(x,t)  this  solution  of  the  difference  equation 
(which  for  noninteger  multiples  t  of  At  is  defined,  for  the  sake  of 

convenience,  as  equal  to  v(x, t*),  t*  =  At[t/At]).  This  v  depends, 
of  course,  on  At  and  Ax. 


A 


Theorem:  Assume  that  as  Ax  and  At  tend  to  zero,  v(x, t) 
converges  boundedly  almost  everywhere  to  some  function  u(x, t) .  Then 
u(x,t)  is  a  weak  solution  of  (1.1)  with  initial  values  cp. 

Proof:  Multiply  Eq.  (1.6'),  satisfied  by  v,  by  any  test  vector 
w,  integrate  with  respect  to  x  and  sum  over  all  values  of  t  which 
are  integer  multiples  of  At.  Apply  summation  by  parts  to  the  left  side, 
and  transform  the  terms  involving  Ag  by  replacing  the  variable  of 
integration  x  by  x  -  Ax/2  and  x  +  Ax/2,  respectively.  We  get  the 
following  identity: 


II 


^,C^t  ,-_At)  -  .w[x,t)  y{X}t)  At 


(1.3') 


-  J w(x,0)  q>( 


x)  dx 


-II 


w(x  +  Ax/2)  -  w(x  -  Ax/2) 


Ax 


g  dx  At; 


g  denotes  g(v^,  . ..,  v^ ^)t  where  v^,  .  ..,  denote  values  of  v 

at  21  points  at  distance  /lx  apart  distributed  symmetrically  around 
x,t.  If  v  tends  boundedly  almost  everywhere  to  u(x,t),  so  do  v^, . .., 

v2i>  and  ‘tllus  S (vi>***>v2/)  ten<is  to  g(u,  . ..,u),  which,  by  the  con¬ 
sistency  requirement  (1.5),  is  f(u).  So  the  limit  of  Cl. 3 * )  is  the 
desired  integral  relation  (1.3). 
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2.  MINIMIZING  THE  TRUNCATION  ERROR 


I 


Let  u(x,t)  be  an  exact  smooth  solution  of  equation  (1.1).  It 

will  then  satisfy  the  difference  Eq.  (1.6')  only  approximately;  the 

deviation  of  the  right  side  from  the  left  side  of  (1.6')  is  called  the 

truncation  error.  It  is  easy  to  see  that  if  the  function  g  satisfies 

2 

the  consistency  condition  (1.5 )>  the  truncation  error  is  at  least  0(A  ). 
In  this  section  we  shall  determine  g  so  that  the  truncation  error  is  of 
as  high  an  order  as  possible .  Specifically  we  shall  consider  the  case 
/  =  1,  and  show  that  g  can  be  chosen  so  that  the  truncation  error  is 
0(A3): 

Expand  u(x,  t+At)  into  a  Taylor  series  up  to  terms  of  second 

order: 


u(x,  t+At)  =  u(x,t)  +  At  ut  +  — +  0 (A3 ) .  (2.1) 

With  the  help  of  the  differential  equation  (1.1)  which  u  is  supposed  to 
satisfy  we  can  express  the  time  derivatives  of  u  as  space  derivatives: 

(2.2) 

10 


ut  “  fx’ 


U. 


=  f  =  f  =  (Au  )  =  (A  u  )  j 

tt  xt  tx  v  t/x  x7x 


what  is  significant  is  that  all  t  derivatives  are  exact  x  derivatives 
and  therefore  can  be  approximated  by  exact  x  differences .  Substitute, 
namely,  (2.2)  into  (2.1);  we  get 


u(x,  tfAt)  =  u  +  (Atf  +  —7"—  A^u  )  +  0(A^) . 

X  X 


(2.3) 


Comparing  (2.3)  with  (1.61),  ve  see  that  the  truncation  error  is  O(A^) 
if  and  only  if 


i  ■  (f+f  A\)x+  0(A2). 


From  this  we  can  easily  determine  the  form  that  g  has  to  take: 

Theorem:  The  truncation  error  in  the  difference  scheme  (1.6*)  is 


0(AJ)  if  and  only  if 


g(a,b)  =  -(a)  -g  -f-^  +  |  A2  •  (b  -  a) 


(2.4) 


plus  terms  which  are  0([a  -  b|  )  for  u  -  v  small. 

Formula  (2.4)  has  a  fairly  intuitive  meaning  and  can  be  derived 

by  differencing  our  differential  equation  (1.1)  as  follows:  replace  time 

and  space  derivatives  by  differences  centered  at  x,  t  +  At/2.  This  means 

that  u,  is  to  be  replaced  by  a  forward  time  difference,  and  f  by 
t  X 


f(x  +  Ax/2,  t  +  At/2)  -  f(x  -  Ax/2,  t  +  At/2)f/Ax. 
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The  value  of  f  at  x  +  Ax/2,  t  +  At/2  is  evaluated  on  the  basis  of 


the  formula 


f (x  ±  Ax/2,  t  +  At/2)  =  f(x  ±  Ax/ 2,  t)  +  ~  ft  +  0(A)2. 

O 

We  express  f  as  Au,  =  A  u  ,  and  approximate  u  by  a  difference 

t)  t)  X  X 

quotient.  The  value  of  f(x  ±  Ax/2,  t)  is  evaluated  as  the  average  of 
f  at  x  and  x  ±  Ax.  The  resulting  formula  is  precisely  (2.4). 

The  quantity  A2  in  (2.4)  shall  be  taken  as  |a2(u)  +  A2(v)j-/2, 
for  sake  of  symmetry  more  than  anything  else;  any  other  choice  for  it 
would  make  a  difference  that  is  quadratic  in  u  -  v . 

Denoting  the  function  (2.4)  by  gQ,  we  can  write  any  permissible 
g  in  the  form 


g  =  gQ  +  |  Q(a>  b)  •  (b  -  a)  ,  (2.5) 

where  Q(a,b)  is  a  matrix  which  vanishes  when  its  two  vector  arguments 
are  equal. 

Substituting  formula  (2.5)  for  g  into  the  difference  Eq.  (1.6'), 

we  get 

2 

u(x,  t+At)  =  u(x, t )  +  XA'f  +  ~  AA^  Au  +  ^AQAu,  (2.6) 

where  A’  denotes  the  operation  |[T (Ax)  -  T(-Ax)],  A  the  operator 


12 


T(Ax/2)  -  T(-Ax/2).  T(s)  denotes  translation  of  the  independent  vari¬ 
able  by  the  amount  s.  We  shall  call  Q  the  artificial  viscosity, 
since  it  occurs  in  the  difference  equations  in  a  way  which  is  similar 
to  the  artificial  viscosity  terms  introduced  by  von  Neumann  and 
Richtmyer  in  [13]  . 

The  difference  Eq.  (2.6)  expresses  the  value  of  u  at  time 
t  +  At  as  a  nonlinear  function  of  u  at  t;  we  shall  denote  this 
function  (operation)  by  N: 

u(t  +  At)  =  N  u(t).  (2.6') 

The  value  of  the  solution  of  the  difference  equation  at  some  later  time 
kAt  is  obtained  from  the  initial  values  by  application  of  the  kth 
power  of  the  operator  N. 

Our  aim  is  to  show  that  the  difference  scheme  (2.6)  is  convergent 
if  the  size  of  \  is  suitably  restricted.  Now  in  the  case  of  linear 
equations  it  is  well  known  and  easy  to  show  (see,  e.g.,  [10])  that  con¬ 
vergence  is  equivalent  to  stability  defined  as  the  uniform  boundedness 
of  all  powers  Nk  of  the  operator  N  within  some  fixed  range  pAt  <  T . 
In  the  nonlinear  case  von  Neumann  has  made  the  reasonable  assumption, 
(see,  e.g.,  [12])  that  the  convergence  of  the  scheme  would  hinge  on  the 
stability  of  the  first  variation  of  the  operator  N.  The  first  varia¬ 
tion  of  N  is  a  linear  difference  operator  with  variable  coefficients; 
von  Neumann  has  conjectured  that  such  an  operation  is  stable  if  and  only 
if  all  the  "localized"  operators  associated  with  it,  i.e.,  the  operators 


obtained  by  replacing  the  variable  coefficients  by  their  value  at  some 
given  point  --  are  stable."^ 

The  stability  of  a  difference  operator  with  constant  coefficients 
is  easily  ascertained  by  making  use  of  the  Fourier  transform.  In  that 
representation  the  application  of  a  difference  operator  with  constant 
coefficients  goes  over  into  multiplication  by  an  "amplification  matrix. 

In  case  the  amplification  matrix  has  n  linearly  independent  eigenvec¬ 
tors,  the  operator  is  stable  if  and  only  if  the  eigenvalues  of  the 
amplification  matrix  do  not  exceed  1  +  O(At)  in  absolute  value. 

On  the  other  hand,  Courant,  Friedrichs  and  Lewy  have  observed  in 
their  classical  paper  [ 2  ]  that  a  necessary  condition  for  the  convergence 
of  a  difference  scheme  is  that  the  rate  of  propagation  of  signals  in  the 
difference  scheme  should  be  at  least  as  large  as  the  true  maximum  signal 
speed,  i.e.. 


Ax 

At 


'max' 


(2.7) 


where  lu  I  is  the  largest  eigenvalue  of  A  at  any  point  within  the 
'^‘max 

relevant  range  of  values . 

Conversely,  we  shall  show: 


1In  the  parabolic  case  the  validity  of  this  conjecture  has  been  estab¬ 
lished  by  Fritz  John  in  his  important  paper  [  4 ]  .  In  the  hyperbolic 
case  a  fragmentary  result  has  be  given  in  Lax  [  8  ] . 


Theorem:  If  condition  (2.7)  is  fulfilled,  the  difference  Eg. 


(2.6)  satisfies  von  Neumann's  condition  of  stability. 

Proof:  The  first  -variation  of  the  operator  appearing  on  the 
right  side  of  (2.6)  is  easily  computed;  its  value  is 

I  +  XAA'  +  |  X2A2A2  +  0(Ax);  (2.8) 

where  A'  and  A  are  defined  as  before,  and  0(Ax)  denotes  an  opera¬ 
tor  bounded  in  norm  by  0(Ax),  provided  that  we  are  perturbing  in  the 
neighborhood  of  a  smoothly  varying  solution,  i.e.,  one  where  neighboring 
values  differ  by  O(Ax).  In  this  case  the  influence  of  the  additional 
artificial  viscosity  term  is  0(Ax);  in  section  3  we  shall  show  how  to 
take  into  account  the  effect  on  stability  of  the  Q  term  in  regions  where 
u  does  vary  rapidly. 

To  "localize"  the  operator  (2.8)  we  merely  replace  the  variable 
matrix  A  by  its  value  at  some  point.  After  making  a  Fourier  trans¬ 
formation,  the  operator  A'  becomes  multiplication  by  i  sin  a,  and 

a  12 

the  operator  A  multiplication  by  2i  sin  — ,  so  that  ^A  becomes 
multiplication  by  cos  a  -  1;  here  a.  is  |Ax,  £  the  dual  variable. 
So  the  amplication  matrix  of  the  operator  (2.8)  is 

I  +  i  sin  a  \k  +  (cos  a  -  1)  X2  A2  +  O(Ax) .  (2.9) 

Denote  the  eigenvalues  of  XA  by  k;  then,  according  to  the  spectral 


mapping  theorem,  the  eigenvalues  v  of  the  amplification  matrix  (2.9) 
are 

v  =  1  +  i  k  sin  a  +  k2 (cos  a  -  1)  +  O(Ax) . 

Since  k  s  \jji  is  real,  the  absolute  value  of  v  is  given  "by 

|v|2  =  (1  -  [1  -  COS  £]  k2)2  +  k2  sin2  £  +  0(Ax) 

=  1  -  2[1  -  cos  £]  k2  +  [1  -  cos  £]2  k4  +  k2 (1  -  cos2  £)  +  0(Ax) 

=  1  -  [1  -  cos  £]2  (k2  -  k4)  +  0(Ax) . 

According  to  our  basic  stability  assumption  (2.7)  the  quantity  k  does 
not  exceed  1  in  absolute  value;  the  above  formula  for  |v|  shows  that 
|v|  is  bounded  by  1  +  0(Ax).  Thus  we  have  demonstrated  that  the 
eigenvalues  of  the  amplification  matrix  do  not  exceed  1  +  0(Ax)  in 
absolute  value,  as  required  in  von  Neumann’s  condition. 
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3 .  ARTIFICIAL  VISCOSITY  FOR  SINGLE  CONSERVATION  LAWS 

In  the  last  section  we  have  determined,  the  function  g(u,v)  up  to 
quadratic  terms  in  u  -  v.  These  undetermined  quadratic  terms  influence 
neither  the  order  of  the  truncation  error  nor  the  stability  of  the  scheme 
at  points  where  the  solution  varies  smoothly.  At  points,  however,  where 
the  solution  varies  rapidly  —  across  a  shock,  that  is  --  it  is  reasonable 
to  expect  that  the  quadratic  terms  have  a  controlling  influence.  In  this 
section  we  shall  show  how  to  choose  the  additional  quadratic  term  Q  in 
our  formula  (2.6)  so  as  to  assure  a  fairly  narrow  shock  width. 

As  observed  in  the  last  section,  Q  enters  the  difference  Eq. 

(2.6)  like  an  artificial  viscosity.  Therefore,  in  order  to  insure  that 
Q  has  a  stabilizing  influence  on  the  difference  equation,  it  is  reason¬ 
able  to  require  that:  a.  Q  is  positive.  Another  property  of  Q,  already 
previously  noted,  is:  b.  Q(a,b)  =  0  when  a  =  b,  while  on  dimensional 
grounds  we  must  require  that:  c.  Q  has  the  dimension  of  A.  We  shall 
discuss  first  the  scalar  case,  i  ,e .,  when  the  number  of  components  of 
u  (and  f)  is  one;  correspondingly  the  matrix  A  =  grad  f  is  one-by-one. 
In  this  case  the  restrictions  a  -  c  show  that  Q  must  be  a  multiple  of 
|A(a)  -  A(b)  | : 
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Q(a,b) 


=  |  |A(a)  -  A(b)|, 


(3.1) 


where  B  is  a  dimensionless  constant  (which  of  course  could  depend  on 
the  value  of  u  and  v) .  With  the  above  choice  for  Q  in  the  differ¬ 
ence  equation  (3.1)  we  shall  now  study  the  shape  of  a  steady  state  shock, 
i.e.,  a  time  independent  solution  of  (2.6)  connecting  states  u(-»)  and 
l4»): 

A'f  +  ^  AA.^  Au  +  ^  AQAu  =  0.  (3*2) 

We  shall  denote  the  values  of  u  at  three  successive  lattice  points  by 
a,  b  and  c.  We  write, 

A*f  =  |[f(c)  -  f(b)]  +  |[f(b)  -  f(a)] 

and  make  the  following  approximation: 

f(c)  -  f(b)  »  A(b,c)(c-b), 
f (b)  -  f (a)  «  A(a,b)(b-a), 

where  A(u,v)  abbreviates  f(u)  +  A(vlt  For  the  simplest  nonlinear 
function  f,  a  quadratic  one,  these  approximations  involve  no  error. 
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Substituting  the  above  changes  into  (3*2 ),  we  get,  after  multi¬ 
plication  by  2: 

|a  +  XA2  +  qJ-  (c-b)  +  |a  -  XA2  -  qJ-  (b-a)  =  0.  (3.3) 

2 

We  make  now  a  further  approximation  by  omitting  the  terms  XA  in  (3 .3), 
under  the  reasonable  assumption  that  the  presence  or  absence  of  this  term 
won’t  alter  much  the  nature  of  the  solution,  since  this  term  has  the  same 
character  as  the  retained  Q  terras,  and  is  expected  to  be  much  smaller 
than  the  Q  term  in  the  important  region  of  rapid  variation,  especially 
if  X.  is  small. 

This  assumption  was  put  to  the  following  numerical  test:  in  the 

2  2 

difference  equation  (2.6),  A  was  replaced  by  WA  ,  where  W(a,b)  is 

a  factor  so  contrived  that  its  value  is  near  1  when  u  and  v  are  close, 

2 

while  its  value  is  near  zero  when  u  and  v  differ  greatly.  The  pres¬ 

ence  or  absence  of  such  a  factor  W  made  hardly  any  difference  in  the 
shape  of  the  steady  state  shock  profile  obtained  experimentally. 

Making  this  change,  and  substituting  our  choice  (3.1)  for  Q  into 
(3-3)y  we  get 


2I.e.,  W  switches  off  the  higher  order  correction  term  in  the  shock 
zone,  where  it  has  no  function  to  fulfill  anyway. 
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|A(b,c)  +  |  |A(c)  -  A(b)|j  (c-b) 
When  B  =  1,  this  equation  can 


+  -|A(a,b)  -  |  |A(b)  -  A(a)|j-  (b-a)  =  0. 
be  written  in  the  following  simple  form: 


Max 


■jA(b),  A(c)j- 


(c-b)  +  Min 


|A(a),  A(b)| 


(b-a)  =  0 . 


(3.3*) 


Denote  by  uL  and  the  two  prescribed  states  at 

are  supposed  to  be  connected  by  a  solution  of  (3*3* )  • 
to  satisfy  the  Rankine-Hugoniot  condition  (1.4)  which, 
shock,  is 


-  oo  and  +oo  which 
These  states  have 
for  a  stationary 


f(uL)  =  f^). 


(3.4) 


In  addition  we  suppose  that  the  entropy  condition  is  satisfied,  i.e*, 
that  the  shock  speed  is  less  than  -A(u^),  the  sound  speed  to  the  left 
of  the  shock,  and  greater  than  -A(u^),  the  sound  speed  to  the  right: 

A^)  <  0  <A(ur).  (3.5) 

It  follows  from  (3*5)  that  there  exists  a  value  u^  such  that 

A(ujjj)  =  o.  (3.6) 
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We  claim  now  that  the  following  lattice  function  is  a  solution  of  (3.3*): 


for  all  negative  lattice  points. 


u(x)  -  <  ^  for  x  =  0, 


(3.7) 


"R 


for  all  positive  lattice  points . 


We  have  to  verify  that  the  difference  equation  (3.3*)  is  satisfied  for 
all  triplets  of  three  successive  lattice  points.  Since  (3.3*)  is  triv¬ 
ially  satisfied  whenever  a  =  h  =  c,  only  three  cases  remain  to  be 
checked: 


I 

II 

III 

a  =  \ 

“m 

b  =  Ul 

C  "  ”m 

"r 

\ 

In  case  I,  it  follows  from  the  first  half  of  inequality  (3*5)  and  from 
(3.6)  that  MaxjA(b),  A(c)j-  =  0;  since  b-a  is  likewise  zero,  (3.3') 
is  satisfied.  We  can  show  similarly  that  (3*3*)  is  satisfied  in  case  II. 
In  case  II,  the  left  side  of  (3*3')  is 


A(V  [UR  •  V  +  A(uL>  Ium  -  V’ 
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Since  A^)  =  °>  this  expression  is  twice 


A(Um>  +  A(Uj,)  r  ,  A<V  .  A(Uj^)  r 

- 2 -  tuR  "  “M1  +  2  ‘  “L1 


We  rewrite  the  two  terms  above  as 


[f(V  -  f(V]  +  [f(uM}  -  f(uL)] 

which  is  correct  if  f  is  a  quadratic  function  and  which  for  general 
nonlinear  f  involves  the  same  kind  of  error  which  we  have  already  com¬ 
mitted.  This  last  expression  is  equal  to  f (u^)  -  f (u^),  a  quantity 
equal  to  zero  by  virtue  of  the  Rankine-Hugoniot  relation  (3*^-)«  This 

verifies  the  difference  equation  (3  •3")  in  case  II* 

We  summarize  our  result  in  this 

Theorem:  Any  two  states  uL  and  Up  which  can  be  connected  by 
a  stationary  shock,  i.e.,  which  satisfy  the  Rankine-Hugoniot  relation 
(3.4)  and  the  entropy  condition  (3*5)*  can  be  connected  by  a  steady  state 
solution  of  our  difference  equation  (3.2).  'This  steady  state  solution  is 
given  approximately  by  formula  (3*7).  which  shows  that  the  transition 
region  is  spread  over  two  mesh  widths. 

The  theorem  above  refers  to  the  case  when  B  is  taken  to  be  1. 

It  is  not  so  easy  to  solve  the  steady  state  difference  equation  for  any 
other  value  of  B,  but  we  conjecture  that  a  steady  state  solution  exists 
for  a  reasonable  range  of  B.  On  the  basis  of  previous  experience  of 
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other  investigators  with  artificial  viscosity  one  would  expect  the  width 
of  the  region  across  which  the  bulk  of  the  transition  takes  place  to 
increase  with  increasing  B.  This  prediction  was  indeed  borne  out  by 
numerical  experiments  recorded  in  the  tables  at  the  end  of  this  paper. 

Likewise,  when  and  u^  are  two  states  which  can  be  connected 

by  a  shock  proceeding  at  some  nonzero  speed,  we  expect  that  the  difference 
equation  (2.6),  with  Q  being  given  by  (3*l)>  possesses  a  steady,  progress¬ 
ing  solution  connecting  these  two  states,  whose  structure  is  similar  to 
the  stationary  solution. 

We  shall  analyse  now  the  stability  of  the  difference  scheme  (2.6), 
with  Q  given  by  (3.1).  We  make  the  assumption  that  the  stability  of 
(2.6)  is  governed  by  the  stability  of  the  following  associated  linear 
difference  operator: 


l  P  p  P  1  2 

I  +  XAA*  +  ^  X  A  A  +  g  , 


(3.8) 


where  A  and  Q  denote  values  of  A  and  Q  at  some  point.  The  ampli¬ 
fication  factor  associated  with  this  operator  is 


1  +  iXA  sin  a  + 


|  jx2  A2  +  Xq|  -jcos  a  -  l|. 


(3.9) 


It  is  not  hard  to  show,  by  an  analysis  parallel  to  that  given  at  the  end 
of  section  2,  that  this  factor  does  not  exceed  1  in  absolute  value  for 
any  frequency  a  if  and  only  if 
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it  follows  from  Eq. 


Denote  "by  |a  the  largest  possible  value  of  |A|j 
(3.1)  for  Q  that  then  Q  will  not  exceed  |bv,  where  V  is  the  var¬ 
iation  of  A .  If  A  does  not  change  sign,  V  <  [i .  Substituting  the 
above  bounds  into  (3. 10),  we  get 

\2  |i2  +  |  B  \  [i  <  1,  (3.11) 

which  is  equivalent  to 

£  (1  +  //16)1/2  -  B/4.  (3.11') 

In  particular,  we  get  for  B  =  1:  — 

<  .  78  • 


This  stability  condition  is  slightly  more  stringent  than  the 
Courant-Friedrichs-Levy  condition  (2  •  7 )  •  It  hears  a  strong  resemblance 
to  the  stability  condition  derived  by  von  Neumann  and  Richtmyer  in  [  13] , 
and  to  the  stability  criterion  derived  later  by  George  White  of  LASL. 

If  in  (3.10)  the  strict  inequality  holds,  then  the  amplification 
factor  (3.9)  is  actually  less  than  1  in  absolute  value  for  all  frequencies 


a  except  a  =  2nn,  n  integer.  This  would  lead  us  to  believe  that  then 
the  steady  state  solutions  are  strongly  stable,  in  the  following  sense: 

Every  solution  of  the  difference  equation  (2.6),  with  Q  given 
hy  (3.1),  whose  initial  values  tend  to  u^  and  u^  as  x  tends  to 
±co,  approaches  a  steady  state  solution  with  increasing  t,  provided 
that  the  stability  condition  (3.11*)  is  satisfied. 

Given  an  arbitrary  initial  state,  such  that  the  corresponding 
exact  solution  of  the  differential  equation  (1.1)  contains  a  number  of 
shocks,  not  necessarily  stationary  or  steady,  progressing,  we  would 
nevertheless  expect  the  corresponding  solution  of  the  difference  equa¬ 
tion  (2.6),  with  the  same  initial  data,  to  bridge  these  shocks  by 
transitions  similar  to  the  stationary  solution  we  have  found  before,  i.e., 
we  expect  the  bulk  of  the  transition  to  be  spread  over  2-3  mesh  widths. 
This  is  reasonable,  since  the  rate  of  variation  of  the  states  at  the  two 
sides  of  the  shock  will  in  general  be  much  slower  than  the  rapid  varia¬ 
tion  within  the  shock  itself. 

Test  calculations  performed  so  far  have  confirmed  this  expectation. 
As  yet  none  of  these  calculations  have  included  very  rapidly  changing 
shocks,  so  at  this  time  we  don’t  know  how  the  present  method  will  handle 
such  a  situation. 
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k.  ARTIFICIAL  VISCOSITY  FOR  SYSTEMS 
OF  CONSERVATION  LAWS 


In  case  of  a  system  of  any  number  of  conservation  laws  the  three 
properties  of  Q  listed  at  the  beginning  of  section  3  no  longer  suffice 
to  determine  Q  within  a  single  dimensionless  constant  but  still  leave 
a  rather  bewildering  variety  of  possibilities. 

Richtmyer  and  von  Neumann  based  their  choice  of  the  artificial 
viscosity  on  a  physical  analogy,  and  the  various  ingenious  modifications 
proposed  and  tested  by  Landshoff,  Harlow  and  Longley  were  likewise  partly 
based  on  physical  analogies.  In  this  section  we  propose  a  form  for  Q 
which  is  dictated  entirely  by  mathematical  analysis . 

Again  we  study  the  shape  of  stationary  solutions  of  (2.6): 

A*f  +  |  x  A  A2  A  +  |  A  Q  A  u  =  0.  (^-1) 


We  denote  again  by  a,  b  and  c  three  consecutive  values  of  u,  and  as 
before  make  the  following  approximation: 


2A‘f  =  f(c)  -  f (a)  =  f(c)  -  f (b)  +  f (b)  -  f(a) 
«  A(b,c) (c-b)  +  A(a,b)(b-a), 


(4.2) 
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Substituting  this 


where  as  before  A(u,v) 
approximation  into  (4.1), 


abbreviate s  ^(u)  g  ^(v)  # 
we  transform  it  to  the  following  form: 


{a+xa2+q} 


(C  -  b)  + 


(b  -  a) 


0, 


(4.3) 


where  A  and  Q  in  the  first  brace  are  to  be  evaluated  between  the 
points  b  and  c,  in  the  second  brace  between  the  points  a  and  b. 

We  wish  to  reduce  the  above  difference  equations  for  a  vector  valued 
function  to  a  scalar  equation.  Such  a  reduction  is  rigorously  possible 
if  all  the  coefficient  matrices  commute;  in  that  case  we  can  get  n 
difference  equations  similar  to  (3.3)  where  the  role  of  A  and  Q  is 
played  by  the  eigenvalues  of  A  and  Q.  This  dictates  the  following 
choice  for  Q:  Q(a,b)  should  be  a  matrix  commuting  with  A(a,b)  whose 
eigenvalues  are  equal  to  the  absolute  values  of  the  differences  of  the 
corresponding  eigenvalue  of  A(u)  and  A(v)  times  dimensionless  factors 
B^,  Bn  of  the  order  of  magnitude  1. 

The  requirement  that  Q  should  commute  with  A  implies,  accord¬ 
ing  to  a  well  known  theorem  of  matrix  theory,  that  Q  is  a  function  of 
A.  This  function  we  can  take  to  be  a  polynomial  of  degree  n  -  1: 


Q 


so  1  +  si  A  + 


(4.4) 


whose  coefficients  gQ,  ...,  gn_1  are  uniquely  determined  by  the  pro¬ 
posed  choice  for  the  eigenvalues  of  Q .  Finding  the  coefficients 
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leads  to  the  Lagrange  Interpolation  problem  which  is 


g0>  •••>  Sn_i 

easily  solved. 

In  the  next  section  we  shall  actually  carry  out  the  determination 
of  Q  in  case  of  the  equations  of  motion  of  a  compressible  fluid  using 
Lagrange  coordinates. 

It  should  be  pointed  out  that  even  if  Q(u,v)  is  chosen  to  commute 
with  A(a,b),  there  is  an  error  involved  in  replacing  (4.3)  "by  scalar 
equations  since  the  matrices  A(a,b)  and  A(b>c)  do  not  commute  in 
general.  Still  we  feel  that  the  above  choice  for  Q  comes  pretty  close 
to  imitating  the  scalar  case . 

We  close  this  section  by  a  further  observation  on  the  shape  of  the 
solutions  of  the  steady  state  equation  (4.1) •  In  analysing  this  equation 
in  the  scalar  case  we  have  simplified  matters  by  dropping  the  term 
which  was  thought  to  be  small  compared  to  the  artificial  viscosity  term 
in  the  region  of  rapid  variation •  On  the  other  hand ^  m  determining  the 
shape  of  the  tail  ends  of  the  transition  curve  the  role  of  the  two  terms 
is  reversed.  We  shall  determine  now  asymptotically  the  shape  of  the 
tail  ends  .  We  write 


a  =  u  +  r  w. 


k+1 

b  =  u  +  r  w. 


k+2 

c  =  u  +  r  w. 
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where  u  is  one  of  the  two  values 


or  Uj  and  w  some  vector 

k 


independent  of  k.  We  substitute  this  into  (4.3),  divide  by  r  and 
drop  all  terms  that  still  contain  the  factor  r  .  In  particular  the 
quadratic  term  Q  drops  out .  We  get 


where  A  denotes  A(u) .  This  equation  has  a  nontrivial  solution  if 
and  only  if  the  matrix  acting  on  w  has  zero  as  eigenvalue.  Now  ac¬ 
cording  to  the  spectral  mapping  theorem,  the  eigenvalues  of  that  matrix 
are 


jd  +  X  d2j-  (r2  -  r)  +  -jd  -  \  d2j-  (r  -  1), 

where  d  stands  for  an  eigenvalue  of  A.  Setting  the  above  expression 
equal  to  zero  gives,  after  eliminating  the  uninteresting  root  r  =  1, 
the  following  value  for  r: 


According  to  the  Courant-Fredrichs-Lewy  condition  (2.7);  |*4i|  1?  80 

that  the  above  expression  is  always  negative .  Furthermore,  at  the  left 
endpoint  u^  we  must  have  |r|  >  1,  at  the  right  end  Up  we  must  have 
| r |  <  1.  This  is  the  case  if  and  only  if  dO^)  is  negative,  dO^)  is 
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positive.  But  according  to  the  entropy  condition  described  at  the  begin¬ 
ning  of  section  1,  two  states  and  u^  which  can  be  connected  by  a 

stationary  shock  always  possess  such  eigenvalues. 

A  similar  analysis  can  be  made  of  the  asymptotic  shape  of  the 
tail  of  steady,  progressing  solutions  of  the  difference  equation  (2.6). 

We  are  led  to  an  equation  which  contains  noninteger  powers  of  r  and 
■which  may  have  complete  solutions. 

All  this  is  not  very  important  as  a  practical  consideration,  but 
it  does  explain  a  curious  phenomenon  observed  already  in  calculations 
performed  with  the  Richtmyer-von  Neumann  method,  and  in  many  calculations 
using  modifications  of  that  method:  that  the  shock  transition  overshoots 
and  approaches  a  constant  value  in  an  oscillatory  fashion.  By  changing 
the  available  parameter,  the  amount  of  overshooting  could  be  diminished 
but  never  completely  eliminated,  nor  could  the  oscillatory  nature  of  the 
approach  to  a  constant  value  be  changed.  The  present  analysis  explains 
this  behaviour  by  the  negative  value  of  r. 

In  contrast,  calculations  performed  by  a  somewhat  crude  earlier 
method  proposed  by  Lax  (see  [6]  and  [7])  produced  shock  transitions 
which  approach  their  final  values  monotonically .  A  similar  asymptotic 
study  of  the  tail  end  of  the  steady  state  solutions  for  these  equations 
disclosed  that  all  the  relevant  values  of  v  are  positive. 
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5.  APPLICATION  TO  HYDRODYNAMICS 


In  section  4  we  have  sketched  a  method  for  constructing  an  arti¬ 
ficial  viscosity  term  for  use  in  a  numerical  scheme  for  arbitrary  systems 
of  conservation  laws .  In  this  section  we  shall  carry  out  the  details  of 
this  construction  in  the  special  case  of  the  equations  of  compressible 
flow  in  Lagrange  coordinates . 

As  dependent  variables  we  shall  use  specific^  volume,  momentum  and 
total  energy,  denoted  by  V,  v  and  E.  The  quantity  v  is  of  course 
velocity,  and  the  total  energy  E  is  the  sum  of  internal  and  kinetic 
energy: 


E  =  e  +  |  v2.  (5.1) 

The  internal  energy  e  is  related  to  pressure  p  and  specific  volume 
V  by  the  equation  of  state 


P  =  P(e,V). 


5I.e.,  per 


unit  mass . 


The  equations  of  conservation  of  mass,  momentum  and  energy  are: 


vt  -  -  v  <5,2) 

Et  -  •  (vp)x- 


Here  t  is  time  and  x  is  Lagrange  mass  variable.  Differentiating 
these  equations  with  respect  to  time,  we  obtain 


Vtt  “  vtx> 
vtt  =  ‘  ptx' 

Et  - 

Using  the  chain  rule,  we  have 

pt  =  pe\+pvV 

According  to  a  well  known  identity  in  thermodynamics 

PPe  -  Pv  =  C, 


(5.3) 


(5» 


(5.5) 
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where  C  is  the  Lagrangian  sound  speed.  We  can  also  write 

et  =  Et"vV  (5*6) 

Using  (5.4),  (5.5),  (5.6)  and  the  original  differential  equations  (5*2) 
to  express  first  t  derivatives  as  x  derivatives,  we  can  rewrite  (5. 5) 
as  follows: 


^tt  “  ^xx’ 


vtt  =  <C  vx>x’ 


(5.7) 


e.  .  =  (pp  +  cr  w  )  . 

tt  '"x  xx 


Recalling  the  identity  valid  for  arbitrary  systems  of  conservation  laws; 


utt  '  <A  ux’x’ 


we  see  from  (5 .7 )  that  in  our  case 


/  - 


A2  u 


C2  v„ 


\  ppx  +  C  wx 


(5.8) 
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Accordingly  we  shall  use  the  approximation 


A2  A  u 


2 

C  A  v 

_2 

pAp+C  v  A  v. 


(5.9) 


We  turn  now  to  the  determination  of  Q,  which,  according  to  the  recipe 
given  in  section  4,  is  a  quadratic  polynomial  in  A(a,b)  whose  eigen¬ 
values  are  constant  multiples  of  the  absolute  values  of  the  differences 
of  the  eigenvalues  of  A (a)  and  A(b).  Now  the  eigenvalues  of  A  are: 
0,  ±  C .  We  claim  that  Q,  is  given  by 


B  |C(a)  -  C (b )  |  .2. 
=  2  C2 


9*  "b 

A  and  C  in  the  above  formula  are  to  be  evaluated  at  ^  •  Clearly, 

Q  as  given  by  the  above  formula  has  the  appropriate  eigenvalues • 

Substituting  this  form  of  Q  into  formulas  (2.4)  (2.5)  for  g  and 
making  use  of  (5.9)>we  have 


I '-*{*%*"}*  \ 

-5+|{fpi+452AT 


(5.10) 
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where  the  barred,  quantities  are  to  be  evaluated  as  averages  between  points 
a  and  b . 

Two  observations  about  this  formula  are  in  order:  first, 

1.  Consider  an  initial  distribution  in  which  v  and  p  are 
constant,  although  V  not  necessarily;  in  fact,  V  may  be  discontinuous . 

For  such  initial  values  the  function  g,  given  by  formula  (5.10),  is  a 
constant;  therefore  the  corresponding  solution  of  the  difference  equation 
(5.2)  is  stationary.  This  agrees,  of  course,  with  the  fact  that  such  an 
initial  configuration  is --in  the  absence  of  diffusion  or  heat  conduction- -in 
hydrodynamic  equilibrium  and  shows  that  our  difference  scheme  introduces 

no  artificial  diffusion  or  heat  conduction. 

More  generally,  we  expect  that  even  in  nonstationary  flows  contact 
discontinuities  are  transmitted  as  sharp  discontinuities . 

2.  Although  our  formula  (5. 10)  for  g  was  derived  on  a  purely 
mathematical  basis,  it  can  be  given  a  more  intuitive  interpretation. 

First  of  all,  as  already  observed  in  section  2,  the  second  order  correc¬ 
tion  terms  can  be  regarded  as  merely  centering  the  values  of  f  properly. 
The  additional  artificial  viscosity  term  can  be  given  the  following 
interpretation:  Define  an  artificial  velocity  and  an  artificial  viscous 
pressure  as  follows: 


mx 


|cj  pv> 


’art  ^2  '"x' 


Mx  1  _  1 

^art  ”  5  l°*l  V 
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where  the  indicated  derivatives  are  to  be  replaced  by  centered  difference 
quotients .  Clearly,  the  effect  of  the  Q  term  in  the  mass  and  momentum 
equations  is  to  augment  the  value  of  velocity  and  pres sure --properly 
centered- -by  the  amounts  indicated  above.  The  effect  of  the  Q  term  on 
the  energy  equation  can  be  described  similarly,  if  we  neglect  the  differ¬ 
ence  between  vp  and  v  p,  and  if  we  neglect  the  product  of  the  artificial 
velocity  and  the  artificial  viscous  pressure. 

The  artificial  viscous  pressure  that  has  cropped  up  in  this  treat¬ 
ment  bears  a  strong  resemblance  to  the  one  introduced  by  Rolf  Landshoff . 

We  present  the  results  of  two  calculations  using  (1.6)  and  (5.10) 
to  obtain  approximate  solutions  of  (5.2).  In  the  first  calculation  we 
initially  had  two  constant  states  separated  by  a  shock.  The  exact  shock 
speed  is  1.  In  Table  I  we  show  the  appearance  of  the  configuration  at  the 
40th  time  cycle,  which,  with  ^  =  .337,  corresponds  to  t  =  15.5.  At 
this  time  the  shock,  which  started  at  x  =  50,  should  be  at  x  =  63*5. 

We  used  2  for  the  parameter  B,  and  the  second  order  correction  terms 
were  not  switched  off. 

In  the  second  calculation  we  again  initially  have  two  constant 
states  separated  by  a  discontinuity.  However,  this  time  the  configuration 
at  40  cycles  is  a  shock  moving  with  speed  1.2b,  a  contact  discontinuity 
at  the  point  of  the  initial  discontinuity,  namely  x  =  50,  and  a  rarefac¬ 
tion  wave.  We  used  ~  =  *337>  B  =  1.  The  results  of  the  calculation  are 

/  \a 

listed  in  Table  II. 
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TABLE  I 

PROGRESSING  SHOCK 


v 


X 

VOLUME 

VELOCITY 

44 

1.000 

1.000 

1.000 

1.000 

46 

1.000 

1.000 

47 

1.000 

1.000 

48 

1.000 

1.000 

49 

1.000 

1.000 

50 

1.001 

1.000 

51 

1.044 

1.000 

52 

1.012 

1.000 

55 

1.005 

1.001 

54 

1.001 

1.001 

55 

1.002 

.998 

56 

1.002 

.998 

57 

.996 

1.006 

58 

•995 

1.007 

59 

1.010 

.986 

60 

1.014 

.980 

61 

.979 

1.052 

62 

.965 

1.054 

65 

1.105 

.864 

64 

1.505 

.411 

65 

1.876 

.085 

66 

1.985 

.011 

67 

1.998 

.001 

68 

2.000 

.000 

69 

2.000 

.000 

70 

2.000 

.000 

ENERGY 

4 .429 

4.429 

4.429 

4.429 

4.429 

4 .428 

4 .434  - - Initial  Position  of  Shock 

4.601 
4.476 
4.445 
4.436 
4 .426 
4.425 
4.442 
4.444 
4.401 
4.388 

4 .496 
4.545 

4 .172  «  Present  Position  of  Shock 

3.584 

2.944 

2 .867 
2.858 

2.857 

2.857 

2.857 
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TABLE  II 


PROGRESSING  SHOCK,  STATIONARY  CONTACT 
DISCONTINUITY,  AND  RAREFACTION  WAVE 


X 

VOLUME 

VELOCITY 

ENERGY 

24 

2 .245 

.698 

20.04 

25 

2.246 

.699 

20  .04 

26 

2.248 

.702 

20  .03 

27 

.2  .253 

.709 

20  .02 

28 

2  .264 

.725 

20.00 

29 

2.284 

.754 

19.95 

50 

2.316 

.800 

19.87 

31 

2.362 

.866 

19.78 

32 

2.422 

.948 

19.66 

33 

2.495 

1.045 

19.53 

34 

2.577 

1.150 

19.40 

55 

2.660 

1.259 

19 .28 

36 

2.756 

1.366 

19.18 

57 

2  .842 

1.463 

19.09 

58 

2.912 

1.541 

19.03 

59 

2.956 

1.589 

19 .00 

40 

2.963 

1.596 

19.00 

41 

2.935 

1 .5  66 

19 .02 

42 

2  .898 

1.525 

19.05 

43 

2.882 

I.508 

19 .06 

44 

2.892 

1.518 

19.05 

45 

2.907 

1.534 

19 .04 

46 

2.909 

1.536 

19.04 

47 

2.903 

1.529 

19.05 

48 

2.902 

1.526 

19.06 

49 

2  .907 

1.528 

19 .08 

50 

2  .887 

1.528 

I8.96 

51 

.825 

I.528 

6.253 

52 

.777 

I.528 

5.956 

55 

.772 

1.528 

5.930 

34 

?774 

1.528 

5-937 

55 

.773 

I.528 

5.932 

56 

.772 

1.528 

5.924 

57 

.771 

1.527 

5.914 

58 

.769 

1.528 

5.906 

59 

.768 

1.527 

5.9OO 

60 

.767 

1.527 

5.895 

61 

.766 

1.528 

5.894 

62 

.765 

1.533 

5.903 

63 

.766 

1.526 

5.889 

64 

.770 

1.519 

5.870 

65 

•749 

1.576 

6.011 

66 

.754 

1.546 

5.983 

67 

1.204 

.850 

4.309 

68 

1.852 

.108 

2.979 

69 

1.990 

.006 

2.863 

70 

2.000 

.000 

2.857 

p  =  3.528 


Rarefaction  Wave 

Exact  Volume  =  2 . 900 


p  =  2  .465 

Initial  Discontinuity 

Exact  Volume  =  .7^7 
Exact  Velocity  =  1.52  8 


Present  Exact 
Position  of  Shock 

p  =  .5714 
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